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ABSTRACT 

The SILCC project (Simulating the Life-Cycle of molecular Clouds) aims at a more 
self-consistent understanding of the interstellar medium (ISM) on small scales and 
its link to galaxy evolution. We present three-dimensional (magneto) hydro dynamic 
simulations of the ISM in a vertically stratified box including self-gravity, an external 
potential due to the stellar component of the galactic disc, and stellar feedback in the 
form of an interstellar radiation field and supernovae (SNe). The cooling of the gas is 
based on a chemical network that follows the abundances of H+, H, H 2 , C+, and CO 
and takes shielding into account consistently. We vary the SN feedback by comparing 
different SN rates, clustering and different positioning, in particular SNe in density 
peaks and at random positions, which has a major impact on the dynamics. Only for 
random SN positions the energy is injected in sufficiently low-density environments to 
reduce energy losses and enhance the effective kinetic coupling of the SNe with the gas. 
This leads to more realistic velocity dispersions (cthi ~ O.Sctsoo-soook ^ 10—20 km s“^, 
~ 0.6(78000-3x105K ^ 20 — 30kms“^), and strong outflows with mass loading 
factors (ratio of outflow to star formation rate) of up to 10 even for solar neighbourhood 
conditions. Clustered SNe abet the onset of outflows compared to individual SNe but 
do not influence the net outflow rate. The outflows do not contain any molecular gas 
and are mainly composed of atomic hydrogen. The bulk of the outflowing mass is dense 
{p ^ 10“^^ — 10“^^ gcm“^) and slow (i; ^ 20 — 40kms“^) but there is a high-velocity 
tail of up to i; ^ 500kms“^ with p ^ 10“^^ — 10“^^gcm“^. 

Key words: hydrodynamics - magnetic fields - methods: numerical - ISM: general 
- ISM: kinematics and dynamics - galaxies: ISM 


1 INTRODUCTION 

Star formation and the resulting stellar feedback, together 
with the thermal, chemical and dynamical evolution of the 
galactic disc, drive the ‘matter cycle’ in the interstellar 
medium (ISM). It describes the way in which gas is cycled 
from a hot, diffuse ionized phase into colder, denser phases, 
and ultimately into molecular gas and back. Star formation 
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occurs within this molecular gas, and the resulting stellar 
feedback both heats the ISM and returns stellar material 
to the ISM. At the same time the cold gas is heated and 
dispersed, while the momentum deposited by outflows and 
winds drives turbulent motions in the gas. The strongest 
feedback is expected to stem from the thermal and dynam¬ 
ical impact of supernovae (SNe), which play a major role in 
the overall dynamics of t he ISM ([McKee fc Ostrik^|1977 


Mac Low Klessen||2004 Klessen Glover ||20 14). 


The effective impact of SN feedback strongly depends 
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on the environment in which the SNe explode ([Taylor 


I950| jSedovI I959| McKee & Ostriker 

I977[ [Cowie, McKee 

& Ostriker|fl98I| [Ostriker & McKee 

19881 

Cioffi, McKee 

& Bertschinger[[I988[ [Slavin & Cox 

I 992 I 

Dwarkadas & 

Gruszko|20I2[ Rogers L Pittard|20I3 

Walch 

& Naab|2015| 

Kim & Ostriker|20I5[[Gatto et al.|20I5[[Iffrig & Hennebelle 


2015 [Martizzi, Faucher-Giguere Qiiataert||2015| |Li et al. 


2015). SNe in dense regions violently interact with their im¬ 


mediate surroundings, potentially destroying molecular gas 
in their vicinity. Efficient cooling can convert most of the en¬ 
ergy into radiation that escapes the SN site and reduces the 
region of dynamical influence. The shells of SNe exploding 
in low-density environments can expand to larger distances 
and reshape the ISM on larger scales. 

Observations of SN remnants in the Milky Way indi¬ 
cate that ~ 20% of the SNe are of type la. Their explosion 
sites follow an exponential distribution centred around the 
midplane with a scale height of ~ 300 pc. The remaining 
SNe are of type II, again with an exponential distribution 
around the midplane but significantly smaller scale height 
of ~ 50 pc (Tammann, Loeffier & Schroeder |I994 ). Stars 
typically form in clusters over a wide range of masses and 
spatial scales [Lada Lada (e.g. 2003). Most of the massive 
stars (> 3/4 of the O stars, de Wit et al.|2004 ) are also clus¬ 
tered, and hence so are most of the type II SNe. However, 
a fraction of the type II SNe are distributed, due to either 


isolated star formation ( 

Kennicutt, Edgar & Hodge I989[ 

McKee & Williams[[I997 

[Clarke & Oey[[2002[ Schilbach & 

Roser|2008[) or runaway stars (Blaauw|I96I Gies & Bolton 

1986 

Gvaramadze & Bomans[2008[[Eldridge, Langer & Tout 

2011 

Perets & Subr|20I2). 


Stellar feedback is also expected to drive galactic foun¬ 
tains, outflows and winds (see, e.g. Chevalier Clegg|I985 


Murray, Quataert & Thompson 2005| ). Observations re¬ 
veal that gas is ejected from the disc at velocities of a 
few 100 km with mass loading factors (ratio of outflowing 
to star-forming gas) of order unity and above (see, e.g. pi I 


et al. 1120 10 


[Martin et al.|2012|[N e 


■ et al.|20I2| [Martin 


et al.||20I3). Simulations of galaxy evolution in a cosmolog¬ 


ical framework support the dynamical impact of fountains 
and outflows and emphasize their importance for metal en¬ 
richment ( Oppenheimer Dave||2008| Oppenheimer et al. 


20I0|[SHnson et al.|20I0||Hopkins, Quataert Murray|20I2| 

Hirschmann et al.||20I3|). 


The properties and the dynamical evolution of the ISM 
have been investigated with numerical simulations of varying 
complexity in terms of the physical processes considered. A 
classical setup to numerically study ISM properties are strat¬ 
ified boxes covering a statistically significant volume of the 
ISM with typical sizes {xxyx ±z) ofx = y~0.5 — I kpc, up 


to z ±20kpc (de Avillez|2000[ de Avillez & Berry|200I[ de 

Avillez & Breitschwerdt |2004[[2005[ Joung & Mac Low 

2006 

Piontek & Ostriker[[2007[ Joung, Mac Low & Bryan[ 

2009 

Koyama & Ostriker[[2009a|b[ Kim, Kim & Ostrikerj 

2011 

Shetty & Ostriker |20 1 2[ 

Hill et al.pOI2[[Creasey, Theuns &[ 

Bower[[20I3[ [Gent et al. 

[20I3a|b[ [Hennebelle & Iffrig[[20I4[ 

Walch et al.|20I5). In the ISM it is important to distinguish 


between the scales of the dynamical evolution of the disc, 
namely the local turbulent and thermal structures, and the 
scales of fountains, outflows, and chimneys that determine 
the large scale gas cycle. For the small scale dynamics in 
the disc, most studies find converged velocity dispersions 


on time-scales of several tens of Myr. The coupling of small 
scale motions to large scale outflows takes an order of magni¬ 
tude longer. Early numerical work on the SN-driven ISM by 
de Avillez & Berry (2001) report the formation of collimated 
chimneys of outflowing gas after ~ 100 Myr. Bothjde Avillezj 
& Breitschwerdt (2004 2005) and Hill et al. ( 20I2| ) empha- 
size that time-scales of a few hundred Myr are needed for 
a fountain cycle to establish a global dynamical equilibrium 
and a converged vertical prohle of the stratified disc. jPio-j 
ntek & Ostriker (2007) use galactic shear driven magneto- 


rotational instability instead of SNe to drive turbulent mo¬ 
tions and establish a dynamical equilibrium in the vertical 
stratification, however, with perceptibly lower velocity dis¬ 
persions in the gas. In their models it takes more than a 
Gyr to establish dynamical equilibrium. jCreasey, Theuns~fcl 
Bower (2013) quantitatively connect the outflow rates to the 


gas content in the disc and the amount of stellar feedback 
which is proportional to the star formation efficiency. They 
find mass loading factors (77 = Moutflow/Afspn) ranging from 
?7 O.I — 4 that have a weak dependence on the star for¬ 
mation rate but significant dependence on the gas fraction 
and surface density of the disc. 

We also note that the Galactic ISM is permeated by 
magnetic fields. On galactic scales the field is expected to 
be generated by the mean-field dynamo with field strengths 
of ~ I — 10 fiG (see e.g. Beck 2001 2009 Grutcher et al. 


2010). The field is dominated by the toroidal component. 
Local tangled components are of similar strengths, however, 
in dense molecular clouds fields up to several mG have been 
measured (see e.g. Grutcher et al.||20I0 Grutcher]|20I2 ). 

Two different aspects of the ISM have not been inves¬ 
tigated in previous studies of stratified box models, namely 
the positioning and clustering of SNe in concert with a more 
sophisticated thermodynamic treatment using a chemical 
network and including shielding and self-shielding of the 
gas. Within the SILGG project (Si mulating the Life-G ycle 
of molecular GloudQ, starting with Walch et al. (2015), we 
aim to advance our understanding of the physical processes 
in the ISM by investigating the impact of differently placed 
SNe as individual explosions as well as clustered explosions. 


Walch et al. (2015) hnd that only simulations with ran¬ 


dom or clustered SN positioning are in agreement with ob¬ 
servations. In those models molecular hydrogen contributes 
40% — 60% to the total mass, whereas the other half of the 
mass is in atomic hydrogen. This paper focuses on the dy¬ 
namical evolution of the gas in the disc. This includes the 
velocity dispersion, the midplane pressures and the resulting 
onset of outflows. 


2 NUMERICAL METHOD AND SIMULATION 
SETUP 

2.1 Numerical Code 

A detailed description of the simulation setup, the numeri¬ 
cal methods and the physical processes covered in this study 
are given in 


Walch et al. (2015). We therefore only pro¬ 


vide a brief description here. The simulations are carried 


^ For movies of the simulations and download of selected simula¬ 
tion data see the SILCC website: www. astro. uni-koeln. de/silcc 
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SILCC: Dynamics of the SN-driven ISM 3 


out with the astrophysical code FLASH in version 4 ( Fryx- 
ell et al.|2000 Dubey et al.|2008 2013), which is an adaptive 
mesh refinement (AMR) code developed at the University 
of Chicago. We solve the magneto-hydrodynamic (MHD) 
equations using the five-wave Bouchut MHD solver HLL5R 


( Bouchut, Klingenberg Waagan|20Q7 2010 Waagan|2009| 


Waagan, Federrath & Klingenberg 2011). We numerically 


solve the following set of equations: 

dp 


^+V-(pv) = 0 (1) 


d{p-v) 

dt 


+ V- 


dE ^ 


T B2\ bb^ 

pvv +(p+_jl-_ 


^ P\ (B-v)B 

E+ — + - V- ^ ' 

Stv p J An 


PV • g Itchem “h ftinj 

<0B 

— - V X (v X B) = 0. 


= PS + qinj 

( 2 ) 

( 3 ) 

( 4 ) 


The total energy density is given by 


E = u-\- 


pv 


B^ 


( 5 ) 


Here, p is the mass density, v the velocity, P — (7 — l)u 
the thermal pressure with u being the internal energy den¬ 
sity and 7 the adiabatic index of the gas, which we set to 
5/3. We do not account for any softening of the equation of 
state in molecular-dominated regions due to the influence of 
the rotational and vibrational degrees of freedom of the hy¬ 
drogen molecule because this effect becomes important only 
when the temperature is high enough to populate the excited 
states. This requires a temperature of at least 100 — 200 K, 
but almost all of the molecular gas we find in our simulations 


is cooler than this (Walch et al. 2015). The magnetic field 


strength is denoted with B. The gravitational acceleration g 
combines self-gravity and the effects of an external potential 
as described below. The total energy density E is coupled 
to the change in internal energy density brought about by 
radiative and chemical heating and cooling, itchem- Finally, 
qinj and itinj describe the mechanical and thermal injection 
rates per unit volume from SNe. 


2.2 Gravity 

For the gravitational forces we take into account self-gravity 
and a background potential due to the stellar component of 
the disc, g = gsg + gext- We include self-gravity of the gas 
by solving the Poisson equation. 


A4> = AivGp, 


( 6 ) 


with a tree based method (Wiinsch et ah, in preparation), 
where 4> is the gravitational potential, G is Newton’s con¬ 
stant and p is the gas density. The external gravitational 
acceleration due to the stellar component in the galactic 
disc is modelled with an isothermal sheet originally pro¬ 


posed by Spitzer (1942), where the distribution function of 
stars is Maxwellian. For our setups we choose a stellar sur¬ 
face density of S* = 3 OM 0 pc“^ comparable to the solar 
neighbourhood (Flynn et al. poo^ |Bo 7^ Rix fc Hoggpol^ 
and a vertical scale height of Zd = 100 pc. We neglect any 
gravitational effects from dark matter because the expected 


local dark matter density (pdm ~ lO~^M 0 pc“^) is more 
than an order of magnitude lower than our average density 
(gas plus stars) in a central volume of (500 pc)^, which is 
(pgas + p*} ~ O.O 5 M 0 pc“^. In addition, the characteris¬ 
tic scale height of the dark matter profile is of the order of 
lOkpc ( Chemin, de Blok Mamon||2011 ). We also neglect 
variations of the stellar potential in the x and y direction 
in order not to impose an initial pattern or a characteristic 
length scale during the dynamical and chemical evolution of 
the gas in the disc. 


2.3 Chemistry and cooling 


We model the chemistry of the ISM using a simplified 
but fast chemical network based on IGlover Mac Fowl 
( 2007a|b ), Glover et al. (2010), Nelson & Langer (1997), and 
Glover Glark ( 2012 ). This network is designed to follow 
the chemical abundances of six species: H^, H, H 2 , G^, GO 
and free electrons; note that for simplicity we assume that 
the abundance of neutral atomic carbon (G) is negligible in 
comparison to that of G^ or GO. 

The thermal evolution of the gas is modelled using the 
detailed atomic and molecular cooling function described 


Glover et al. (2010), Glover & Glark ( 2012 ) and Walch 


et al. (2015). This cooling function makes use of the infor¬ 


mation on the chemical composition of the gas provided by 
our chemical network, and hence is significantly more accu¬ 
rate than the simple analytical functions used in a number 
of previous studies ( Walch et al.||2011 Micic et ar]|2013| . 

We assume a uniform interstellar radiation field (ISRF), 
which we scale linearly with the number of massive stars and 


thus the SN rate (see Sec. 2.4). For our fiducial mod els with 
a SN rate of 15Myr“^ we use a value of Gq = 1.7 (jDraine 


1978) in units of the Habing field, Go ( Habing||1968 ). The 


ISRF is attenuated in dense regions due to gas and dust 
shielding based on the column density using the TreeGol al¬ 
gorithm ( Glark, Glover Klessen||2012 ). Full details of the 
implementation of this chemical network in FLASH 4 and 
the way in which it is coupled to our TreeGol-based treat¬ 
ment of molecular self-shielding and dust shielding can be 


found in Walch et al. (2015) and Wiinsch et al. (in prepara¬ 
tion). 

In the simulations presented in this paper, we assume 
solar metallicity, with fixed elemental abundances of carbon, 
oxygen and silicon given by xc = 1.41 x 10 ~^, xp — 3.16 x 
10“^ and xsi = 1.5 x 10“^, respectively (Sembach et al. 


2000). Initially, we start with all of the carbon in the form 


of G^ and all of the oxygen in the form of O. Silicon is 
present in the form of Si^, and as it does not participate 
in the chemical network, is assumed to remain in this form 
throughout the simulation. 


2.4 SN energy input 

Feedback is provided by SNe, which is injected as thermal 
energy or as momentum input depending on the resolution 
and the density of the injection region. If the radius at the 
end of the Sedov-Taylor expansion phase (i.e. the radius at 
which the SN expansion changes from the non-radiative ex¬ 
pansion phase to the radiative snowplough phase) is resolved 
with at least 4 cells, we inject an energy of 10^^ erg per SN 
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4 Girichidis et al. 


explosion as thermal energy and let the code convert ther¬ 
mal into kinetic energy self-consistently. If the density in 
the volume under consideration is too high to resolve the 
Sedov-Taylor expansion phase with pcrit = 3 x 10“^^ gcm“^ 
in our simulations, injecting thermal energy would lead to 
over-cooling because the resulting net temperature in the 
injection region would be below ~ 10® K, where the cool¬ 
ing rates are very high. The SNe would then have a neg¬ 


ligible dynamical and thermal effect (see, e.g. Anninos & 
Norman||1994 Stinson et al.||2006 [Creasey et al. |2011| ). In- 

stead, we inject the expected net momentum the gas would 
be exposed to if the substructures of the SN explosion were 
resolved. The physical details we employ are described in 


Blondin et al. (1998). For a detailed investigation of the SN 


positioning and energy injection methods in a very similar 
chemical environment we refer to Gatto et al. (2015). 


The SNe are injected at a constant rate from the very 
beginning of the simulation. From the Kennicutt-Schmidt 
relation ( [Schmidt] [1959 Kennicutt 1998| hereafter simply 
called KS relation) we derive a star formation rate (SFR) 
for our simulation box, which in turn can be converted into 
an SN rate using a fixed initial stellar mass function. We 
assume the formation of one massive star per 100 Mq of 
stars. This yields a SN rate of 15 Myr“^. We account for the 


variations and uncertainties in the KS relation (Bigiel et al. 
2008||S^ruba et al.|201 11 [Leroy et al.|2013||Shetty, Kelly 

Bigiel|2013[ [Shetty et al.j2014[ ) by varying the SN rate by a 

factor of three above and below the standard value. 

The positions of the SNe are chosen in different ways. 
We place them either in the density peaks (peak driving), 
at random positions in the xy-plane and a Gaussian dis¬ 
tribution in z with a standard deviation of 50 pc (random 
driving), or alternating between the two modes with half of 
the SNe exploding in density peaks and half at random posi¬ 
tions (mixed driving). The random positions are determined 
beforehand and stored in a SN sequence to ensure the same 
positions in all simulations which use the random driving 
mode and the KS SN rate. In addition, we perform three 
runs with clustered SNe. 

For the runs with clustered SNe we define different ver¬ 
tical scale heights for SNe of type la and II, simila r to [dej 
Avillez & Breitschwerdt (2004) or Joung & Mac Low ( 2006[ ). 
We split the total number of SNe into a fraction of 20% 
type la SNe with a scale height of 325 pc and 80% SNe of 


type II with a scale height of 50pc (Miller & Scalo 1979 
Heiles|19^ Tammann, Loeffler Scliroeder|1994 ). Among 

the type II SNe, 3/5 explode in a clustered environment, the 
rest of the type II SNe explode as individual SNe at random 
positions. For the clustered SNe we assume a total cluster 
lifetime of 40 Myr, draw the number o f SNe per cluster, iV, 
from a power-law distribution oc N~‘^ (Kennicutt, Edgar & 


Hodge[[T^ [McKee Williams[[T^ (Glarke Oey[[2002 ) 

with a minimum of A" = 7 and a maximum of A = 40 SNe 
per cluster, and spread the A SNe with equal temporal ex¬ 
plosion intervals dt — 40Myr/A over the cluster’s lifetime 
( Joung Mac Low|2006[ ). Each cluster has a fixed position 
over its lifetime. When the next SN in the stratified box 
is due we determine the kind of explosion (type la, type II 
in cluster, type II individual) according to their fractional 
abundance, and if it is a cluster type II SN we find the one of 
the predefined clusters, whose next SN is closest to the simu¬ 
lation time. After all SNe of a cluster exploded the cluster is 


removed from the list. New clusters with new random posi¬ 
tions are created (or better activated) if the current SN that 
is due according to the KS rate does not fit in the explosion 
rate of one of the existing active clusters. 

In addition we perform a run with the KS SN rate but 
with all SNe clustered as type II with a scale height of 50 pc 
as used in the case of random driving. This run allows us 
to distinguish between effects that come from SN clustering 
and effects that are due to type la SNe. 

The justification of SNe at random positions - and 
therefore most likely in low-density environments (see Sec.[^ 
- is based on several aspects. Eirst, there is a non-negligible 
fraction of runaway OB stars (e.g. de Wit et al.[[2005 ) that 
can travel a few hundred parsecs before they explode as a 
SN. Eor our simulation box with an area of 0.5 x 0.5 kpc^ 
in the xy-plane a random positioning seems reasonable (see 
discussion in Li et al. 2015). In addition, recent work by 


Bressert et al. (2012) and Oey et al. (2013) suggests that 


some of OB stars is born in the field. Another important 
aspect is the numerical resolution. A SN that explodes in a 
low-density environment in our simulations would need to 
travel at least a few cells away from the imaginary high- 
density birth place of the progenitor star. At a resolution 
of 4 pc this distance can easily be ~ 40 pc. Higher resolu¬ 
tion would result in more resolved structures including cav¬ 
ities and voids next to dense filaments and cores, which are 
smoothed to 4 pc in our setups. Statistically, the same star 
would therefore need to travel much less in order to reach 
low densities. Besides the problems of the travel distance 
of massive stars we also expect them to create low-density 
regions around them early on via protostellar outflows and 
later via radiation and stellar winds. This does not justify 
a random position, but it justifies the accompanying low- 
density environment (see also Gatto et ah, in preparation). 


2.5 Initial conditions 

The initial density profile of the ISM in our model is uniform 
in the xy-plane and follows a Gaussian distribution in the 
vertical direction. 


p{z) = 1 ^ 


( 7 ) 


with zo = 30 pc and po = 9 x We apply a 

density floor of pmin = to mimic the hot at- 

mosphere around the disc. The resulting column density of 
E = 10 Mq pc“^ is similar to the Galactic value at the so¬ 
lar radius, E 13M0pc“^ (Elynn et al. 2006). The gas 


is initially at rest. We deliberately refrain from setting ini¬ 
tial turbulent motions and density perturbations in order 
not to introduce specific scales or dynamical patterns. In¬ 
stead, the SN driving together with the chemical evolution 
and the dynamics due to gravitational forces is allowed to 
form structures and motions self-consistently. The gas in the 
disc is initially atomic at a temperature of 4500 K. The hot 
regions above and below the plane consist of ionized gas 
with a temperature of T = 4 x 10^ K. We set the molec¬ 
ular gas fraction to zero everywhere in the box to follow 
the formation history of H 2 . The disc is set up in thermal 


pressure equilibrium (Wang et al. 2010). The presence of 


self-gravity and external potential thus result in a disc that 
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Table 1. Simulation parameters 


Quantity 

Symbol 

Value 

box size 

X X y X z 

0.5 X 0.5 X ±5kpc 

effective resolution 

Ax 

4 pc 

gas surface density 

E 

10 Mq pc~^ 

gas disc std dev 

ZQ 

30 pc 

central gas density 

PO 

9 X 10-24 gcm-3 

stellar surface density 

S* 

30 Mq pc-2 

stellar disc scale height 

Zd 

100 pc 

std dev. SN type la 

Zsma 

325 pc 

std dev. SN type II 

^^SNII 

50 pc 


is initially not in hydrostatic equilibrium. However, test sim¬ 
ulations where we switch off cooling and chemical evolution 
show that gravitational effects only marginally change the 
density profile while reaching hydrostatic equilibrium. The 
effects of cooling are much more severe in driving the disc 
out of equilibrium. 

In the magnetic run we assume that the field is initially 
oriented in the x direction. The field strength scales with 
the density, 

B{z) = (8) 

V Po 

with the field strength at z = 0 being Bq = 3 /xG. No tangled 
component is introduced. 


2.6 Simulation overview 

The overall numerical parameters are listed in Table ^ Our 
simulated box spans 500 pc in x and y direction and extends 
to ±5kpc in z direction. The effective resolution is set to 
Ax = 4 pc, corresponding to 128 x 128 cells in the xy-plane 
for 1^1 < 2kpc. For |z| > 2kpc we reduce the resolution to 
Ax = 8 pc. Table [T6| gives a list of all simulations discussed 
in the paper. The simulation name is composed of the surface 
density of the disc (S = 10 Mqpc“^ ^ SIO), the SN rate 
where lowSN, KS and highSN refer to 1/3, one and three 
times the KS SN rate. The SN driving mode is indicated 
by rand (individual randomly placed SNe), peak (individual 
SNe in density peaks), mix (mixed mode with half the SNe 
in random places and half in density peaks) and clus (clus¬ 
tered SNe with random positions of the clusters). In run 
clus2 all SNe are of type II and clustered without single in¬ 
dividual SN explosions. The suffix nsg denotes a comparison 
run without self-gravity. The simulation including magnetic 
fields with an ordered field strength in x direction of 3 jaG 
in the midplane has mag attached to the simulation name. 


3 MORPHOLOGICAL EVOLUTION 

In Fig. we show a comparison of the face-on (integrated 
along the z axis) column density after lOMyr of evolution 
for runs with a KS star formation rate. At this early phase 
SlO-KS-peak (left panel) shows clear signatures of the peak 
driving algorithm. The first SN creates an overdense shell 
where the subsequent SNe are placed automatically. This 
coherent clustered SN input sweeps up a lot of the gas 
and it takes about ~ 40Myr to completely mix the gas 


in the disc. For random SN driving (SlO-KS-rand, second 
panel in Fig. the early morphology is very different. The 
gas develops a filamentary structure with column densities 
exceeding O.lgcm”^ (densities exceeding 


g cm 


m 

-2 


(10“^^gcm“^) in Voids’. This early structure is created 
mainly by SN explosions. Self-gravity has little effect at this 
evolutionary stage (see SlO-KS-rand-nsg). With clustered 
driving (SlO-KS-clus) a similar filamentary morphology de¬ 
velops, however with larger structures caused by the larger 
number of spatially coherent SN explosions. The presence 
of the magnetic field delays the formation of dense struc¬ 
tures. After lOMyr the magnetic field run (SlO-KS-clus- 
mag3, right panel in Fig.[^ shows less well developed struc¬ 
tures at lower column densities. 

In Fig. we show the edge-on and face-on morphology 
of the simulations after t = 30Myr for simulations SIO-KS- 
peak (top) and SlO-KS-rand (bottom). We plot (from left 
to right) a density and temperature cut through the posi¬ 
tion of the maximum density in the box followed by pro¬ 
jections of the total gas density, H^, H, H 2 , and CO. For 
peak driving (SlO-KS-peak) the disc is more compact. The 
molecular component is dilute and most of the gas has inter¬ 
mediate and low temperatures (< 10^ K). The accumulation 
of gas induced by the peak driving algorithm is clearly vis¬ 
ible in the face-on plots. The structure of the disc does not 
change substantially during the evolution and no outflows 
are launched. 

Already after 30 Myr the disc is significantly thicker (in 
particular for neutral and ionized hydrogen) with random 
supernova driving (SlO-KS-rand) and the densities and tem¬ 
peratures indicate a noticeably larger volume filling fraction 
of hot (> 10^ K) gas near the disc plane (for a quantita¬ 
tive analysis see Walch et al. 2015). The same qualitative 
behaviour can be seen for simulations with random driving 
and low and high star formation rates. The expanding neu¬ 
tral and ionized components indicate the onset of an ISM 
outflow as discussed in more detail in Sec. |5.4| The H 2 and 
CO columns show a number of dense clumps of molecular 
gas in the disc midplane which by this time already con¬ 
tributes ~ 40 per cent to the total gas content (only ~ 7 per 
cent for SlO-KS-peak, see Walch et al. poT^ . 

The differences between random (SlO-KS-rand) and 
clustered (SlO-KS-clus) driving are highlighted in Fig. at 
a later stage of the simulation at t = 100 Myr. Clustered SN 
explosions have mainly two effects. They drive gas flows to 
larger heights above the disc (described in detail in Sec. 5.4) 
and generate more hot gas with a larger volume filling factor. 
Local fluctuations in the gas distribution combined with the 
random - but statistically centred around z = 0 - position¬ 
ing of the clustered SNe allow for an asymmetric evolution 
with respect to the z = 0 plane. Dense regions start to form 
above and below the plane. Over time the SNe, whose ver¬ 
tical distribution remains centred around z = 0, can thus 
drive gas asymmetrically out of the plane, which explains 
the offset position of the massive cloud. In the subsequent 
analysis we thus define the midplane as the position of the 
densest gas for all simulations except SlO-KS-peak and SlO- 
KS-rand-nsg, whose densest peaks fluctuate perceptibly. In 
the disc region the gas is more concentrated in clumps with 
a larger fraction of H 2 in the run with clustered SNe com¬ 
pared to one with individual random SNe. At t = 100 Myr 
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Table 2. List of simulations 


no. 

simulation name 

SN rate 
(Myr-l) 

/rand 

/clus 

/type II 

B 

(mG) 

self¬ 

gravity 

1 

SlO-lowSN-rand 

5 

1.0 

0.00 

1.0 

- 

yes 

2 

SlO-lowSN-peak 

5 

0.0 

0.00 

1.0 

- 

yes 

3 

SlO-lowSN-mix 

5 

0.5 

0.00 

1.0 

- 

yes 

4 

SlO-KS-rand-nsg 

15 

1.0 

0.00 

1.0 

- 

no 

5 

SlO-KS-rand 

15 

1.0 

0.00 

1.0 

- 

yes 

6 

SlO-KS-peak 

15 

0.0 

0.00 

1.0 

- 

yes 

7 

SlO-KS-mix 

15 

0.5 

0.00 

1.0 

- 

yes 

8 

SlO-KS-clus 

15 

- 

0.48 

0.8 

- 

yes 

9 

S10-KS-clus2 

15 

- 

1.00 

1.0 

- 

yes 

10 

S10-KS-clus-mag3 

15 

- 

0.48 

0.8 

3 

yes 

11 

SlO-highSN-rand 

45 

1.0 

0.00 

1.0 

- 

yes 

12 

SlO-highSN-peak 

45 

0.0 

0.00 

1.0 

- 

yes 

13 

SlO-highSN-mix 

45 

0.5 

0.00 

1.0 

- 

yes 


Column 3 gives the supernova rate per Myr. The SNe are distributed with a certain 
fraction of random locations, /rand? which is given in column 4: /rand = 1-0 corre¬ 
sponds to purely random SN driving, whereas /rand = 0-0 corresponds to pure peak 
SN driving. Column 5 and 6 give the fraction of clustered SNe and the fraction of SNe 
of type II, respectively. Column 7 contains the magnetic field strength and column 8 
indicates whether self-gravity is switched on or off. We use the same notation for the 
simulations as |Walch et al.| ( [2QI5| >. 


200 
100 
' 0 
-100 
-200 


SlO-KS-peak 


SlO-KS-rand 


- 200-100 0 100 200 
X (pc) 



S 10-KS-rand-nsg 


SlO-KS-clus 


S 10-KS-clus-mag3 


S 


/ '7 


M . , 7 

7 / 


A _ 


X (pc) 


X (pc) 
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Figure 1. Early evolution of the column density for different SN driving modes at t = 10 Myr. Shown are from left to right SIO-KS- 
peak, SlO-KS-rand, SlO-KS-rand-nsg, SlO-KS-clus and S10-KS-clus-mag3. The run with peak driving creates a superbubble due to the 
subsequent explosions of SNe in swept-up over-densities of the first SN. The forming structures in the clustered SN run are larger than in 
the case of individual random positions of the SNe because of multiple explosions at the same position. Self-gravity does not play a role 
at this stage of the simulation. The presence of magnetic fields causes dense structures to form later because of the stabilising magnetic 
pressure. 


the small molecular clouds that have formed initially in the 
filamentary ISM have merged into a few very massive molec¬ 
ular complexes. A time sequence of this process is shown 
in Fig.j^for SlO-KS-rand. After stable molecular structures 
have formed they are pushed towards each other in the pres¬ 
ence of turbulent motions at scales of > 10^ pc. At smaller 
distances they gravitationally attract each other and merge 
into larger complexes. The merging clouds are not disrupted 
because most of the random SNe and clusters are located 
within the volume-filling hot phase. The morphological evo¬ 
lution indicates that the simulations do not reach an equi¬ 
librium configuration in the morphology (see Sec. for a 
discussion). 

The different morphological evolution of the ISM com¬ 


bined with the SN driving mode results in different envi¬ 
ronmental densities, in which the SNe explode. In Fig. 
we present the distribution of the environmental density, in 
which we inject the SNe. We average over time excluding the 
initial phase of the simulation (t < 20 Myr). The vertical line 
indicates the threshold value for the thermal/kinetic injec¬ 
tion mode. SNe at lower densities are injected as thermal 
energy. For higher densities we switch to a mechanical injec¬ 
tion scheme (for details see Section 2.4). For peak driving 
all SNe turn out to be mechanical. Mixed SN positions show 
a broad bimodal distribution with environmental densities 
ranging from 10~^®gcm“^ to 10~^^gcm“^. The environ¬ 
mental densities are even higher than in the peak driving 
run because only half of the total number of SNe explodes 
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Figure 2. Overview of the simulations SlO-KS-peak (top) and SlO-KS-rand (bottom) ai t = 30Myr. We plot (from left to right) a 
density and temperature cut through the position of the maximum density in the box followed by projections of the total gas density, 
H+, H, H 2 , and CO. The y position of the cuts in the top panel is indicated by the white line in the bottom panel. Similarly, the 2 : 
position at which the cut in the lower panel is taken is indicated by the white line in the upper one. For peak driving (SlO-KS-peak) 
the disc is very compact. The molecular component is dilute and most of the gas has intermediate and low temperatures (< 10"^ K). For 
random SNe (SlO-KS-rand) the midplane contains more warm and hot gas. The density structures are clumpy and more molecular gas 
has formed. 
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SlO-KS-clus 



10-1 10 “^ 10 “^ 10-1 ]_ o - 


Figure 3. Overview of the simulations SlO-KS-rand (top) and SlO-KS-clus (bottom) at t = 100 Myr. We plot (from left to right) a 
density and temperature cut through the position of the maximum density in the box followed by projections of the total gas density, H+, 
H, H 2 , and CO. The y position of the cuts in the top panel is indicated by the white line in the bottom panel. Similarly, the 2 ; position 
at which the cut in the lower panel is taken is indicated by the white line in the upper one. Initial filaments and small clouds merged to 
eventually form a few dense, massive GMCs. Simulation SlO-KS-rand shows a dense column of mainly atomic gas above and below the 
plane. Temporal statistical asymmetries in the SN positions cause the formation of a CMC at \z\ 100pc in simulation SlO-KS-clus. 

Strong outflows that do not contain molecular gas have developed in both simulations. 
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t = 10 Myr t = 40 Myr 



t = 70 Myr 



t = 100 Myr 



10-5 


10-'‘ 10-5 JO-5 10-5 

Column density (g cm-^) 


Figure 4. Evolution of the total gas column density at t = 10,40, 70,100 Myr for simulation SlO-KS-rand. Over time the morphology 
changes from thin filaments and over multiple diffuse clouds that form at the intersections of filaments. The individual low-mass clouds 
both accumulate more gas and merge with other clumps to form a few very massive molecular clouds at the end of the simulation. 



Figure 5. Distribution of the environmental density of the SNe 
for all simulations with KS SN rate averaged over the simulation 
time for t ^ 20 Myr. Below the critical density, Pcrit? we inject 
thermal energy, while above it we use mechanical injection to 
avoid over-cooling. In run SlO-KS-peak, all of the SNe explode 
in regions denser than pcrit with a peak of the distribution at 
10“22gcm“3, gQ q]\ Qf g]'^g injected as mechanical 
injection scheme. Random and clustered models are dominated by 
thermal explosions in the density range 10“^^ — gcm“^. 


in density peaks, which reduces the net counteracting force 
against gravitational attraction in dense regions. All ran¬ 
dom and clustered models are dominated by thermal energy 
injection in densities from 10 “^^ — 10 “^^gcm“^. 


4 DYNAMICAL EVOLUTION 

To quantify the dynamical evolution of the gas we compute 
its one-dimensional velocity dispersion, as, accounting for 
the thermal cr^,therm and turbulent as,turh contributions: 

(^s,turb “1“ therm (9) 

Here, s indicates the ”gas phase” under consideration. In 
our analysis, we consider two ways to partition the ISM 


into different phases. First, we use the traditional separa¬ 
tion into four temperature regimes, (Ti > 3 x 10^ K, T 2 G 
[8000K; 3 X 10^ K], T 3 e [300K; 8000K], T 4 < 300K), 
which are usually taken to represent the hot ionized medium 
(HIM), the warm neutral medium (WNM), thermally un¬ 
stable gas between the WNM and the cold neutral medium 
(CNM), and finally the CNM itself, which is not separable 
from the even colder molecular gas in this approach. Our 
second method for identifying different ISM phases is more 
physically motivated and closer to observations. It is based 
on our ability to follow the chemical evolution of the gas 
(Sec. |2.3] ). We estimate the flux in the Hi 21 cm line and the 
Hq: line and use these as tracers of the atomic and ionized 
gas, respectively. Ideally, we would use a similar approach 
for the molecular component. However, the resolution in our 
current set of simulations is not sufficient to allow us to make 
accurate predictions for the CO emissivity, and we use the 
following approach for H 2 instead. 

We calculate the mean turbulent velocity dispersion as 


C^s,fc,turb — 


V Ms 


1/2 


( 10 ) 


where k is the spatial direction (x, y, z) and i is the index 
of the computational cell. The velocity Vk,i is the velocity 
of cell i and = (Ei is the mean velocity in 

direction k. When computing the velocity dispersion of H 2 , 
or of the different temperature components, we use simple 
mass weighting. In this case, ms,i represents the mass of the 
relevant phase in cell i, and Ms represents the total mass 
in that phase in the entire computational volume. Alter¬ 
natively, when computing the velocity dispersion of Hi or 
Ha, we weight their contributions proportional to their esti¬ 
mated flux, as explained in more detail below. In this case, 
ms,i represents the flux from cell i and Ms the total flux. 
The mean turbulent one-dimensional velocity dispersion is 
then 


C^s,turb — 


1 ^ 

t 2 
o '5's.fc. 


\ 1/2 
turb I 


( 11 ) 
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Figure 6. Velocity dispersions for all simulations with the KS SN rate. The left column depicts the total velocity dispersion (turbulent 
+ thermal), the middle and the right column present the turbulent and thermal contribution. The top row shows mass-weighted values 
of (7 h2 5 the middle and bottom rows show the observationally motivated estimates for (Jhi and fJHa- The dynamics in the dense gas 
component, crH2, is small for the simulation without self-gravity and for SlO-KS-peak. The evolution of cthi and differs noticeably between 
random and clustered driving. For H 2 and Ha the turbulent component clearly dominates. In the case of Hi, the turbulent component 
only dominates for clustered driving. In all other runs both contributions crturb and crtherm are similar. 


The thermal contribution to the velocity dispersion (thermal 
broadening) is given by 


therm 


Ws 


1/2 


( 12 ) 


with the thermal velocity i;therm,z = the 

temperature, Ti, the mean molecular weight in cell i, and 
the mass of a hydrogen atom, mu- 

The separation into temperature regimes is straightfor¬ 
ward. To estimate the velocity dispersion for the 21 cm emis¬ 
sion from atomic hydrogen, we assume the emission to be 
optically thin. This is a reasonable assumption for all but 


the highest column density lines-of-sight ( 

Heiles & Troland 

2003 Fukui et al.||2014 Motte et al.||2014 

Bihr et al.||2015 ) 


and allows us to determine crHi,turb directly from the atomic 
hydrogen number density provided by the chemical network. 
The Ha intensity is not directly proportional to the mass in 
H+ as the emission decreases with temperature. To com¬ 
pute the Ha fluxes we follow the same procedure as outlined 


Gatto et al. (2015) accounting for collisional excitation 


and radiative recombination of ionized hydrogen (Dong & 
DrainepOll Draine|[2011| Kim et ^^[2013 ). We note that 


although most emission is expected from gas around 10^ K 
there can be contributions from radiative recombination at 
lower temperatures. 

In Fig. we show the time evolution of the velocity 
dispersions crH 2 , ctri, and aua (from top to bottom) for the 
seven simulations with a KS star formation rate. The total 
values (turbulent plus thermal) are on the left, turbulent 
and thermal contributions (see Eq. are given in the mid¬ 
dle and right panels, respectively. We do not distinguish be¬ 
tween intra- and inter-cloud dispersion. Molecular hydrogen 
(top panels of Fig.[^ starts forming from very cold gas only 
after ~ 10 Myr (see |Walch et ah 2015) and the turbulent 
component of the velocity dispersions slowly rise from ~ 2 
to ~ 6kms“^ for all runs except the two with peak driving 
(green line) and random driving without self-gravity (yellow 
line), whose values do not exceed 2.5kms“^ until the end of 
the simulation. The effect of self-gravity is therefore clearly 
visible in the molecular gas kinematics. As the H 2 gas is very 
cold (T < 300 K), the total dispersion is dominated by tur¬ 
bulent motions with negligible contributions from thermal 
broadening. 

The inferred velocity dispersions of neutral hydrogen 
(Hi, middle panels of Fig. slowly rise from initial val- 
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shows the smallest values for both the turbulent as well as 
the thermal components. The strong effect of self-gravity 
seen in the H 2 dynamics is almost absent (yellow line) as Hi 
traces warm and more diffuse gas. For all runs with individ¬ 
ual SN explosions (no clustering) the thermal broadening is 
comparable to the turbulent a, whereas for clustered SNe 
the turbulent contribution clearly dominates. 

Observable Ha dispersions (bottom panels of Fig. are 
roughly constant over time within the scatter with values of 
~ 20—30 km . The mixed and random driving simulations 
show slightly lower values at the end (~ 20 — 25kms“^), 
whereas the clustered driving runs reach 30 —40 km s“^ with 
the highest numbers for the simulation with only clustered 
SNe, S 10 -KS-clus 2 . A noticeable exception is again SIO-KS- 
peak, in which basically all SNe are embedded in dense re¬ 
gions and little hot gas is accelerated. The fluctuations are 
stronger because for this hotter gas phase we start to see the 
direct impact of individual supernova explosions. For Ha dis¬ 
persions the turbulent motions now dominate over the ther¬ 
mal broadening (ana, therm — 12kms“^ for gas around 10^ K 
as shown in the bottom panels of Fig. see also Fig. |^. 
Similar to Hi, there are no strong effects from self-gravity. 
Magnetic fields (run S10-KS-clus-mag3) do not seem to per¬ 
ceptibly change the dynamical evolution of the systems, in¬ 
dependent of the tracer. 

Previous studies of the evolution of stratified discs did 
not include detailed models for the chemical evolution of the 
gas. They typically only focused on gas in different tempera¬ 
ture regimes (e.g. de Avillez Breitschwerdt||2005| ). We use 
what we consider our most realistic model, i.e. the run with 
magnetic field and clustered driving (S10-KS-clus-mag3), to 
investigate which temperature regime is probed best by the 
tracers we have just discussed. In the upper panel of Fig. 
we collect the results from S10-KS-clus-mag3 (red lines in 
Fig. §. In the middle panel we show the evolution of the 
gas dispersions in the five different temperature regimes. 
The hot ionized gas (Ti > 3 x 10^ K) will be emitting in X- 
rays (not investigated here) and reaches a peak dispersion of 
> 200kms“^. The dynamics of the warm ionized gas (crT 2 ) 
is well traced by aua. The same applies to the warm neutral 
gas in the temperature regime T 3 and the Hi measurements 
(see the yellow (Ha/T 2 ) and green (H 1 /T 3 ) lines in the bot¬ 
tom panel of Fig.|^. The velocity disperison of the cold gas 
(T 4 ) correlates well with the values for H 2 , 0 - 112 /(7t4 = 0.4. 
Similarly, the ratio for Hi/Ha shows remarkably small tem¬ 
poral scatter for t > 50Myr with crni/crHa = 0.6. 

Our temperature-based velocity dispersions (middle 
panel of Fig.[7| a gree well with those reported by de Avillez 


& Breitschwerdt (2005), Joung, Mac Low & Bryan ~p009| ). 


Koyama & Ostriker (2009a), Kim, Kim & Ostriker (2011), or 
Shetty & Ostriker ( 2012 ), who all investigated the SN-driven 


ISM. The temperature-binned velocity dispersions found by 


Piontek & Ostriker (2007) are slightly lower than ours. How¬ 


ever, we note that these authors examined the magneto- 
rotational instability as a driver of turbulence rather than 
SNe. 


Our simulations with random and clustered driving 
with self-gravity are also reasonable models for the ob¬ 
served ISM on spatial scales ranging between a few 10 
to 100 parsecs. The dispersion in molecular hydrogen (di¬ 
rectly determined from H 2 in the simulations) agrees well 




Figure 7. Time evolution of the velocity dispersion for simulation 
S10-KS-clus-mag3 for the chemical species and observational esti¬ 
mates (top), split up in temperature regimes (middle) and various 
ratios (bottom). The dynamics obtained from chemical species 
corresponds very well with the dynamics obtained from temper¬ 
ature cuts, indicated by the same colour in the top and middle 
panel. All ratios in the bottom panel show very little temporal 
variation. 


with observed (in CO) molecular gas complexes or giant 
molecular clouds (GMCs) in M31 {a = 6.5 ± 1.2kms“^, 
Sheth et al. 12008]) and GMCs in ot her star formin g galaxies 


(1 — lOkms ^ ,|Bolatto et al.|2008 ). More recently. 


Donovan 

Meyer et al. ( 2Q13| ) reported values of around 8 krns“^ for 


spatially resolved observations of GMCs. To a similar level, 
our molecular dispersions agree with other GMC measure¬ 


ments in M31 (Rosolowsky 2007), the LMC (Fukui et al. 
2008), and the SMC (Mizuno et al. 200 1 | see also the re- 


view by Fukui fc Kawamura||2010 ). lanjamasimanana et al. 


( 2012 ) use data from the THINGS survey (Walter et al. 


2008) not only to estimate dispersions for the cold neutral 
medium {a = 6.5 ± 1.5kms“^) but also for the warm neu- 
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tral medium. Their value of a = 16.8 di 4.3kms“^ agrees 
well with our Hi dispersious. 

lu our simulatious am domiuates over am throughout 
the simulatiou time with crHi/crH 2 ~ 3 — 4. Assumiug that 
<7co = crH 2 this ratio is somewhat larger thau observed val¬ 
ues by |Caldu-Primo et al. ( |2013| ), where most systems show 
crm/crco ~ 1 — 2 . 

The estimates of our velocity dispersious for are 
iu agreemeut with observatious. Greeu et al. (2014) observe 
galaxies with resolved kiuematics fiudiug values iu the rauge 
of ana ~ 15 — 30kms“^ for Milky Way like systems. 


5 VERTICAL STRUCTURE AND OUTFLOWS 
5.1 Vertical density structure 

As the vertical deusity structure is complex, we do uot at¬ 
tempt to ht it with a simple fuuctioual form, but iustead 
determiue the distauce from the midplaue euclosiug 60% 
{Hqo%) aud 90% {Hgo%) of the mass. For au expoueutial 
profile p{z) = po exp(—|z|//i) the total mass iu the rauge 
[—h,h] is about 60%, which motivates the choice of Hqq%. 
As almost all deuse molecular gas is coufiued to the disc 
midplaue the value of Hqq% is couuected to the deuse struc¬ 
tures aud the molecular gas, whereas i/ 90 % also iucludes the 
euvelope of the disc. 

Fig.[^shows the time evolutiou of Hqo% (left) aud Hqo% 
(right) for all simulatious with KS SN rate. Iu both ruus with 
iudividual raudom SNe (SlO-KS-raud aud SlO-KS-raud-usg) 
the iuuer part of the disc (1^60%) quickly expauds. If self¬ 
gravity is switched off the expausiou starts uoticeably earlier 
aud leads to a thicker disc over the eutire simulatiou time. 
Clustered driviug models cau efficieutly compress gas iuto 
deuse filameuts aud clouds, which results iu smaller values 
for Hqq% for the first half of the simulatiou time. Ruus SIO- 
KS-clus aud S10-KS-clus-mag3 eveutually drive euough gas 
out of the midplaue to reach Hqq% ~ 400 pc, i.e. more thau 
40% of the total mass is fiually driveu out of the disc. Iu 
the case of ouly clustered SNe (S 10 -KS-clus 2 ) the efficieut 
compressiou due to the cohereut SNe results iu more mas¬ 
sive agglomerates of deuse gas close to the midplaue, which 
keeps the values for the 60% mass limit below ~ 30 pc over 
the eutire simulatiou time. Mixed SN driviug euds up with 
similarly low uumbers for Hqq% but with a differeut evolu¬ 
tiou. The combiuatiou of half the SNe explodiug iu peaks 
aud more volume filliug raudom uuclustered SNe first cause 
the disc to expaud. The compact gas iu the midplaue ouly 
forms later wheu iudividual deuse clouds merge aud the frac- 
tiou of SNe iu deusity peaks cauuot preveut the clouds from 
mergiug further. Simulatiou SlO-KS-peak with a compact 
disc of mostly diffuse gas shows a very similar evolutiou for 
both Hqo% aud Hqo%. 

Iu the evolutiou of the euvelope of the disc {Hgo%) we 
uotice a clear differeuce betweeu iudividual raudom aud clus¬ 
tered SNe driviug. Ou the oue haud clustered SNe result iu 
more gas close to the midplaue, ou the other haud they cau 
more efficieutly drive au outflow, which causes the euvelope 
to expaud much faster. Duriug the first 50Myr the mixed 
driviug model forms au exteuded diffuse disc, which does uot 
allow the SNe do push gas to larger heights. Iu the secoud 
half of the simulatiou time outflows are eveutually lauuched 
aud the euvelope cau expaud. 


A more detailed view of the disc structure is showu iu 
Fig. where we plot the vertical distributiou of the mass 
fractious as a fuuctiou of time for simulatious SIO-KS-raud- 
usg, SIO-KS-raud, SIO-KS-clus, aud SI0-KS-clus-mag3. We 
also overplot the heights for 60% (solid hue) aud 90% (dot¬ 
ted hue) euclosed mass. The ruus with raudom driviug (up¬ 
per pauels) show smoothly distributed gas with uo stroug 
features iu the profiles aud au overall similar evolutiou of 
Hqo% aud Hgo%. The two ruus with clustered driviug (lower 
pauels) show strouger deusity coutrasts iu the profile aud a 
much larger differeuce betweeu the evolutiou of the expaud- 
iug euvelope {Hgo%) aud the exteut of the deuse part of the 
disc {Heo%). 

Compariug the vertical structure to other models is dif¬ 
ficult because we cauuot ruu the simulatious for loug euough 
to establish a full fouutaiu cycle aud a dyuamical equilibrium 
ou large scales over time. Outflows reshape the profiles uutil 
the eud of the simulatiou. However, H 2 aud CO form aud 
remaiu close to the midplaue of the disc iu all our calcula- 
tious, which allows us to compare those distributious to ob¬ 
servatious. We therefore compare H 2 scale heights with the 
heights of 60% euclosed mass iu H 2 , i^H 2 , 60 %? ciud similarly 
for 60% of euclosed mass iu CO, i^co, 60 %- In most simula¬ 
tious H 2 is very couceutrated with i^H 2 , 60 % < 30 pc except 
for simulatious SIO-KS-raud-usg aud SIO-KS-peak for which 
^H 2 , 60 % slowly iucreases to ~ 100 pc at t = 100 Myr. How¬ 
ever, for those two models the mass tractiou of H 2 is below 
5% aud the molecular gas is very diffuse. We hud CO al¬ 
ways embedded iu the deusest parts of the H 2 clouds aud 
so i^co, 60 % ^ ^H 2 , 60 %- Rcccut work by Lauger, Piueda & 


Velusamy (2014) fiuds a scale height based ou au expoueutial 


fit of 47 pc for CO iu the Milky Way. Observed scale heights 
iu other galaxies cau be as small as 40 pc for CO, but cau 
also reach values of up to 200 pc ( Caldu-Primo et al.||20I3 


Yim et ffi]|20I4 ). A possible reasou for why our molecular 


scale heights are ou the low side of the observed rauge might 
simply be resolutiou. With our curreut resolutiou of 4 pc, we 
are primarily seusitive to the largest GMCs, which predom- 
iuautly form iu the midplaue, aud we might miss the small, 
trausluceut high-latitude clouds that we see iu the Milky 
Way. ludeed, observatious of molecular clouds iu the Milky 
Way show that the scale height of massive GMCs is cousid- 


erably less thau that of smaller molecular clouds (Stark & 
Lee||2005 ), cousisteut with this explauatiou. Also, we would 


like to stress auother limitatiou of our simulatiou setup. Our 
box represeuts ouly a small fractiou of the galactic disc aud 
is uot exposed to large scale galactic dyuamics. The gas iu 
the box evolves iu isolatiou aud is uot iuflueuced by warps 
iu the disc, iucliued accretiou flows, large scale fouutaiu ef¬ 
fects, aud mergers. All of these effects are likely to act as au 
additioual source of motiou that might easily iucrease the 
vertical scale height of deuse gas. Iu additiou, other physical 
processes like CR pressure aud radiative feedback are likely 
to euhauce the dyuamics of the disc aud iucrease the scale 
height of the gas. 

The vertical distributiou of the fractiou of atomic hydro- 
geu as a fuuctiou of time is showu iu Fig. for simulatious 
SIO-KS-raud-usg (top left), SIO-KS-raud (top right), SIO- 
KS-clus (bottom left) aud SIO-KS-clus (bottom right). The 
rest of the mass iu the visible area is maiuly H^ because 
almost all the H 2 is coufiued withiu |z| ~ 100pc. Over¬ 
plotted are the heights for 60% {Hqo%) aud 90% {Hgg%) of 
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Figure 8. Vertical heights of 60% left) and 90% right) enclosed mass over time. SN positions in density peaks result in 

compact discs with small values for Random and clustered SNe lead to a continuous expansion of the envelope {Hqq%) without 

any sign of a slowdown or turnover. Clustered driving models push the gas to noticeably higher altitudes. The height of the innermost 
60% of the gas differs between random and clustered driving modes. Without self-gravity the inner part of the disc starts expanding 
after only 10 Myr. Self-gravity delays the expansion of the innermost 60% of the mass. In the clustered runs the gas is pushed into dense 
structures more efficiently resulting in small values of Hqq% until a significant fraction of the gas is expelled in outflows. For S10-KS-clus2 
more than 60% of the mass are confined in one CMC at the end of the simulation. 



Figure 9. Time evolution of the mass fraction as a function of height for simulations SlO-KS-rand-nsg, SlO-KS-rand, SlO-KS-clus, and 
S10-KS-clus-mag3. Overplotted are the vertical heights of 60% (solid line) and 90% (dotted line) enclosed mass. The simulations with 
random SN positions have relatively massive and smooth vertical gas distributions, in particular at later times. The clustered SN models 
including 20% of type la SNe show low-density regions between the Hqq% and line. 


the total enclosed gas mass. We notice a difference in the 
chemical structure between the random and clustered driv¬ 
ing models. The former runs have smooth atomic hydrogen 
dominated gas distributions with a smoothly decreasing H 
fraction above the Hqq% line. The latter simulations indi¬ 
cate a more complicated composition along the z direction 
mainly due to the clustered SNe and the more powerful hot 


gas chimneys that are created in those simulations. The ad¬ 
ditional type la SNe only show a minor effect, which we 
discuss in more detail in Section [5.41 


Vertical density profiles averaged over 5 Myr of evolu¬ 
tion are shown in Fig. 11 at t = 100Myr (right). The cen¬ 
tral density of the simulation without self-gravity is more 
than two orders of magnitude smaller than in all runs in- 
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t (Myr) 


t (Myr) 


Figure 10. H fraction over time for SlO-KS-rand-nsg (top left), SlO-KS-rand (top right), SlO-KS-clus (bottom left) and SlO-KS-clus- 
mag3 (bottom right). Overplotted are the heights of 60% {Hqq%) and 90% ( 1 ^ 90 %) enclosed mass. As H 2 only forms at very low heights, 
which can hardly be seen in this linear plot, the remaining mass of hydrogen is in the form of H+. The outflows are overall dominated 
by atomic hydrogen. Individual random SNe lead to very high ratios of atomic to ionized hydrogen. In models with clustered SNe the 
composition locally drops to yield roughly equal proportions of H and H+. 


eluding self-gravity. The peak driving model shows very flat 
density profiles up to z ~ 200 pc, where the profile drops 
steeply. The simulations with only clustered SNe can effi¬ 
ciently drive outflows which is reflected in the lower densities 
above z ~ 100 pc. 


5.2 Midplane pressure 

Gas motions are governed by thermal (Ptherm), kinetic 
(^kin = cind magnetic (Pmag = P^/Stt) pressure, 

which are shown in cuts through the midplane for the 
non-magnetic runs in Fig. at t = 30 Myr (top) and 
t = 100 Myr (bottom). At t = 30 Myr the values span five or¬ 
ders of magnitude, so the gas is clearly not in pressure equi¬ 
librium. For all simulations the kinetic pressure dominates 
over the thermal pressure, in particular in dense regions. In 
the subsequent evolution larger and larger structures form 
through merging of filaments and clouds, which is also visi¬ 
ble in the pressure distribution at t = 100 Myr. 

In Fig. we show a comparison of midplane-pressures 
for run S10-KS-clus-mag3. As Pmag is small compared to 
the other two components, we multiply the magnetic pres¬ 
sure by a factor of 10 to better visualise the results. The 
flux-freezing of the ideal MHD approximation leads to a 
correlation of the magnetic pressure with the density. At 
t = 100 Myr the magnetic pressure contribution is confined 
to the dense structures in the midplane. 

A more quantitative analysis is given in Fig. |14[ where 
we show the mean pressure in a thin slice averaged over the 


midplane. The left panels show the values at t = 30 Myr for 
all simulations, the right panels are at t = 100 Myr. In the 
top row we use the values at the midplane, in the bottom 
plots we average the numbers in the z range ±100 pc. The 
black lines are the values directly obtained from the simula¬ 
tion box. The grey area indicates the midplane pressures in¬ 
ferred from observations by Blitz & Rosolowsky (2006). For 
this comparison we have to compute the ’’interstellar pres¬ 
sure”, Pext, in the same way as Blitz & Rosolowsky (2006), 
cf. their equation (1), 


(2G)'/'Egag 



(13) 


with G being Newton’s constant. Eg the surface density of 
the gas, cTg the total velocity dispersion of the gas, and p* 
and pg the stellar and gas densities. For the stellar density, 
p*, we use the same values as for the external potential (see 
Sec. 1^. The values for Pext are shown with the coloured 
lines. 

The numbers support the visual impression that the ki¬ 
netic pressure dominates over the thermal one in almost all 
runs. In the simulations SlO-highSN-mix and SlO-highSN- 
rand the disc is already blown apart which explains the 
low kinetic and high thermal pressures. The total pressures 
taken directly from the midplane in the simulations are at 
the upper end of the observed values at t = 30 Myr. At 
the end of the simulation at t = 100 Myr almost all numbers 
clearly exceed the observed range. Computing Pext using the 
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Figure 12. Thermal and turbulent midplane pressures at t = 30Myr (top) and t = 100Myr (bottom). Overall the kinetic pressure 
dominates over the thermal one. At t = 30 Myr the pressure contrasts are large and reflect small-scale structures. Over time the small 
scale pattern in the pressure distribution evolves into large scale structures along with the coalescence of gas into fewer massive clouds 
towards the end of the simulation. 


gas density and the velocity dispersion yields values which 
are still high but in agreement with the observations. 


Since the observations by Blitz & Rosolowsky (2006) are 
only resolved down to ~ 100 pc it seems reasonable to also 
compare their pressures to the simulation values based on 
averaged quantities from ±100 pc. The central high values 


of the kinetic pressure are then attenuated by a factor of 
up to an order of magnitude. Using the averaged quantities 
from ±100 pc to compute Pext shows a similar agreement 
of the total pressures taken from the simulations with the 
observed estimates. 
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Figure 14. Midplane pressures for all simulations at different times. The upper plots show the pressures in the midplane at t = 30 Myr 
(left) and at t = 100Myr (right). The lower plots show averaged pressures in the 2 ;-range [—100pc; 100pc]. The black lines denote values 
directly taken from the simulation data, split into kinetic and thermal components. The coloured lines are the ’’interstellar pressure” 
computed using Eq. ( |13| ) ( [Blitz fc Rosolowsky|200^ Eq. (1)). The grey area indicates the range of observed midplane pressures by |Bhtz| 
I&: Rosolowsky||2006| >. In almost all cases the kinetic pressure dominates over the thermal one. Using the midplane simulation data yields 
higher values than the estimated ’’interstellar pressure”. Averaging the region within ±100pc reduces the relative difference below a 
factor of two. 


5.3 Vertical pressure profiles 

In Fig.^jwe plot from top to bottom vertical profiles of the 
thermal (Ptherm) and the total pressure, 

Ptot = {pv^} + (Ptherm) + ^ , (14) 

OTT 

as well as their ratio at t = 30 Myr (left) and t = 100 Myr 
(right). For t < 50 Myr the thermal pressure shows moderate 
variations and local fluctuations along the z coordinate. It is 
not surprising that the peak driving run has the lowest val¬ 
ues of Ptherm bccausc the SNe are mostly injected as momen¬ 
tum rather than as thermal energy. The thermal pressure in 
the midplane ranges from Ptherm/fe = 2 — 20 x 10^ cm“^ K. 
At the end of the simulation the thermal pressure at high 


altitudes decreases by an order of magnitude due to cooling. 
The central values do not change significantly over time ex¬ 
cept for run SlO-KS-mix, which forms a dense GMC towards 
the end of the simulation, in which energy is deposited by the 
fraction of SNe placed in density peaks. In contrast the total 
pressure changes noticeably. At t = 30 Myr Ptot drops by a 
factor of 10 from the centre to high altitudes {z > 0.5kpc). 
The ratio of kinetic to thermal pressure (Pkin/Ptherm = 
with the Mach number M) below reveals that the inner¬ 
most ~ 200 — 1000 pc are dominated by kinetic contribu¬ 
tions. Above that height the thermal pressure dominates 
and the flow becomes sub-sonic. Up to that evolutionary 
stage the difference between the individual runs is not dra¬ 
matic (except for the height at which the ratio drops). At 
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Figure 15. Vertical profiles of the thermal (top) and total pressure (middle) as well as the ratio of kinetic to thermal pressure (bottom) 
at t = 30Myr (left) and t = 100Myr (right). The thermal pressure shows moderate variations along the 2 : coordinate at an early stage 
of the simulation. After 100 Myr the variations increase up to two orders of magnitude. On the other hand, the total pressure shows 
an outward gradient early on which becomes stronger over time. At t = 100 Myr the central pressure differs noticeably. Runs SlO-KS- 
peak and SlO-KS-rand-nsg show similar values which are an order of magnitude lower compared to the other runs. At t = 30 Myr the 
ratio of kinetic to thermal pressure indicates that the kinetic pressure is ~ 5 times higher than the thermal one in the disc. Above 
2 : 200 — 1000 pc the thermal component dominates and the flow becomes sub-sonic. 


© 0000 RAS, MNRAS 000 , 000-000 






































18 Girichidis et al. 



Figure 11. Vertical density profiles at t = 100 Myr. The mid¬ 
plane region indicates very similar results with p(z = 0) ~ 
10“23gcm“3 for all except two simulations. The one with peak 
SN driving shows a very compact smooth disc with a flat den¬ 
sity profile up to 200 pc. The central density is about an order 
of magnitude below the values of the other runs. The second no¬ 
ticeable exception is the random driving run without self-gravity 
with more than two orders of magnitude lower central densities. 
The efficiency of the SN clustering in shaping the environment is 
visible at heights above 2 : ~ 100 pc where the density in SIO-KS- 
clus2 is more than an order of magnitude lower compared to the 
runs with a lower clustering fraction. 


t = 100 Myr the total pressure in the midplane increased 
by a factor of a few up to one order of magnitude except 
for run SlO-KS-rand-nsg. This inner part of the disc up to 
z ~ 50 pc is dominated by kinetic pressure and coincides 
with the z positions of the majority of the SNe. The height 
from 50 — 500 pc is at pressure ratios of order unity. The 
drop in total pressure at z > 1 kpc marks the front of the 
outflowing gas. 


5.4 Outflows and mass loading 


Because of the short simulation times we cannot discuss 
galactic outflows and fountain effects in a steady state (see 
de Avillez Breitschwerdt|2Q04 Hill et al.|2012 ) . However, 


we can investigate the initial onset of an outflow as a func¬ 
tion of driving mode and SN rate. We measure outflow ac¬ 
tivity as the flux of mass at 2 : = 0.5 kpc and at 2 ; = 1 kpc and 
relate this quantity to the star formation rate. This gives the 
mass loading factor 


ri{z) = 


M{z) 

Msfk 


(15) 


Outflowing and inflowing gas is traced separately instead 
of only tracing the net flow. We would like to emphasize 
that inflow measured at height 2 : is only gas that has been 
pushed out before and now falls back towards the midplane 
of the disc. We have no external accretion flow onto the disc, 
i.e. from outside the simulation box. In addition, we analyse 
the composition of the outflow by splitting up the outflowing 
mass into the individual hydrogen components. 

The previous analysis demonstrates that the peak driv¬ 
ing model has a compact disc with no or very little mass 
being launched out of the disc. In this section we thus dis¬ 
card this model. In Fig. we summarize the outflow prop¬ 
erties over time. The plots on the left are the values at 




Pressure (cm ^ K) 


Figure 13. Midplane distribution of the thermal (top), kinetic 
(middle) and magnetic pressure (top) for simulation SlO-KS-clus- 
mag3 at t = 30Myr (left) and t = 100Myr (right). The magnetic 
pressure is overall weaker and shown magnified by a factor of 
10. The ideal MHD approximation correlates the density struc¬ 
ture with the field strength (flux freezing). The patches of high 
magnetic pressure thus trace the density field. 


2 ; = ±0.5 kpc, the plots on the right at 2 ; = ±1.0 kpc. The 
top row shows the outflow (solid lines) and inflow (dotted 
lines). The chemical composition is shown as the ratio of 
atomic hydrogen mass over the mass in H^ in the second 
row. None of the outflows carries any H 2 or CO; they are 
composed solely of H and H^. This result is possibly af¬ 
fected by resolution. With a minimum cell size of 4 pc we 
are unable to resolve typical substructures in dense regions. 
Higher resolution allows for the formation of denser regions 
that can form molecular gas more easily. In addition, stable 
dense structures might also resist the dissolution by hot fast 
gas more effectively. The third and fourth row present the 
density and the mass-weighted outflow velocities, separately 
for atomic hydrogen (solid lines) and ionized hydrogen (dot¬ 
ted lines). 

We note that our simulations establish strong outflows 
after a short time of ~ 10 — 30 Myr. As we start our calcula¬ 
tions with a simple isothermal disc without any substructure 
and velocity fluctuations, we need to wait for the SNe to in¬ 
sert enough energy and momentum before measurable out- 
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Figure 16. Time evolution of the outflow properties at 0.5kpc (left) and 1 kpc (right). The upper plots show the outflow (lines) and 
inflow rates (dots) which are generally lower. A coherent outflow is only launched with random or clustered driving and only after 
t ~ 10 — 30 Myr (0.5 kpc) or t ~ 30 — 50 Myr (1.0kpc), respectively. Comparing the star formation rate (SFR, see Sec. |2.4| > to the outflow 
rate yields mass loading factors of up to 10 at the end of the simulations. The panel below depicts the composition of the outflow. We 
note that the outflow does not contain any H 2 but mainly consists of atomic hydrogen. For clustered driving at 2 ; = 0.5 kpc the fraction 
of H is lower in favour of H+. The third row shoes the outflow density, which is remarkably similar for all runs in atomic hydrogen. The 
outflow velocity in the bottom panels is a factor of a few lower for atomic hydrogen compared to ionized hydrogen. 


© 0000 RAS, MNRAS 000 , 000-000 































20 Girichidis et al. 


flows are launched. Clustered SNe start launching outflows 
earlier than individual SNe. The inflow is smaller over the 
entire simulation time such that we have net outflow with 
a mass loading factor of up to 10. Only for random driving 
does the inflow increase towards the end with inflow rates 
larger than the star formation rate. At a height oi z — 1 kpc 
we notice a similar behaviour with a delay of ~ 20 Myr cor¬ 
responding to a net outflow velocity of 25kms“^ 

The ratio of H to in Fig. 


16 


indicates that the com¬ 
position is dominated by atomic hydrogen for all runs except 
S10-KS-clus2 at the end of the simulation. Strong SN clus¬ 
tering generally lowers the amount of atomic hydrogen. The 
fact that S 10 -KS-clus 2 without any SNe of type la shows the 
lowest ratio suggests that the SN positioning/clustering in 
the disc has a stronger impact on the composition of the out¬ 
flow than the inclusion of individual type la SNe with a six 
times larger distribution in z. At z — 1 kpc the correlation 
with the driving mode is less pronounced. 

We compute the mean outflow density assuming that 
within one computational cell the chemical mass fractions 
equal the corresponding volume fractions. At z = 0.5 kpc 
there is only little scatter in the density of atomic hydrogen 
with p ~ 2 — 10 X 10“^^gcm“^ for all simulations, slowly 
rising over time. The values at 1 kpc are very stable with 
similar densities. For ionized hydrogen we measure about 
one to two orders of magnitude lower numbers with a larger 
temporal scatter. 

The panel at the bottom depicts the mean outflow ve¬ 
locities. At z = 0.5 kpc there is a difference between the 
mass-weighted H velocities of about 10kms“^ and the 
velocities which are three times as high at the end of the 
simulation. S10-KS-clus2 shows the largest H velocities of 
~ 50kms“^ at the end of the simulation. At z — 1 kpc the 
scatter in the velocities is very high with a less pronounced 
distinction between H and due to turbulent instabilities 
and mixing. This behaviour is illustrated in Fig. |17| where we 
show the density with velocity vectors (left) and the frac¬ 
tion of the gas (right) in cuts through the box for simulation 
SlO-KS-rand nt t — 100 Myr. We find high velocities of up 
to several hundred kms“^ in the low-density gas that starts 
as a collimated chimney and mixes with the high-density, 
low-velocity gas at z > 0.5 — 1 kpc. 

In Fig.[^we illustrate how the density and the compo¬ 
sition of the outflow vary with the vertical velocity. Plotted 
are two-dimensional histograms for Vz—p (top) and Vz — fii+ 
(bottom) for the outflowing gas for simulation SlO-KS-clus- 
mag3 at t = 100 Myr. The colour-coded outflow rate indi¬ 
cates an almost bimodal distribution, with the bulk of the 
outflow moving at 10 — lOOkrns”^ at a mean density of 
p ~ 10“^^ — 10 “^^gcm“^. Although the high-velocity gas 
{vz > lOOkrns”^) is at low densities it contributes signifi¬ 
cantly to the total outflow rate. 

Given the outflow velocities and densities we can esti¬ 
mate the amount of gas that can escape from the disc for 
our choice of the external gravitational potential neglect¬ 
ing further outward acceleration. In order to escape from 
the Milky Way the gas at the solar radi us must have ve - 
locities above Uesc(lOkpc) ~ 500kms“^ (Xue et al. 2008). 
= 500kms“^ 


For Ue, 


the fraction of escaping mass is only 


2 — 8 x 10 “ ,so basically no gas is expected to escape from 
Milky Way gravitational potential based on the z velocities 
we measure at the end of the simulation. The fraction of 
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Figure 17. Cuts through the box of the gas density with outflow 
velocity vectors (left) and ionization fraction (right) for simu¬ 
lation SlO-KS-rand at t = 100 Myr. At \z\ = 0.5 kpc the out¬ 
flow follows the structure of a collimated low-density chimney, 
at 1 2 ;I = 1 kpc the flow converts into more turbulent motions. 
The low-density ionized gas reaches velocities of a few hundred 
kms“^ whereas the dense gas moves at velocities of the order of 
10 kms“^. 


the volume with outward velocities above 500kms“^ is sig¬ 
nificantly larger and ranges from 1 — 30% based on a total 
volume of 0.5 x 0.5 x 4kpc^ excluding the regions of unper¬ 
turbed pristine gas far above and below the disc. 

Determining the mass loading factor, 77 , from observa¬ 
tions is very difficult and involves many assumptions about 
the composition of the gas and the structure of outflows 
(see discussion on clumpy versus smooth outflows in the ap¬ 


pendix of Martin et al. 2013). Observed outflows indicate 
mass loading factors of order unity ( Chen et al .|2010 


Martin 


et al.| 20 l 2 Newman et al.| 2012 | ). Simulations of galaxies also 
agree with those numbers (see, e.g. Oppenheimer & Dave 


2008 Hopkins, Quataert Murray| 2012 ). At z = 0.5 kpc we 
reach values of up to 77 = 10 which are higher than the ob¬ 
servational estimates. However, we only simulate the onset 
of fountains and winds for time-scales shorter than needed 
for typical galactic fountains. Gas that eventually falls back 
towards the midplane can perceptibly decrease the net mass 
loading factor. A conclusive statement concerning the con¬ 
nection between local mass loading factors at z = 0.5 — 1 kpc 
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Figure 18. Two-dimensional histograms for Vz — p (top) and 
Vz — /h+ (bottom) for the outflowing gas for simulation SIO-KS- 
clus-mag3 at t = 100 Myr. Colour-coded is the outflow rate. There 
is almost a bimodal distribution in the outflow rate (upper plot). 
The bulk of the outflow has velocities of a few tens of kms“^ 
at a mean density of p ~ lO^^'^gcm”^. We note that the high- 
velocity gas {vz > 100kms“^) is fully ionized and mostly at 
very low densities. However, this high-velocity tail contributes 
significantly to the total outflow rate. Below Vz ^ 100kms“^ the 
fraction of ionized gas shows a wide spread. 


and escaping gas {z ^ lOkpc) cannot be drawn based on 
our simulations. 


Creasey, Theuns & Bower (2013) investigate the mass 


loading factor in detail as a function of disc surface density 
and the ratio of gas to dark matter. Their surface densities 
are overall higher than our value of E = 10 Mq pc“^, which 
does not allow for a direct comparison. In addition, they 
impose a temperature floor at 10^ K, which might measur¬ 
ably change the structures in the midplane where the SNe 
explode and thus the outflow properties. However, extrapo¬ 
lating their curves of the mass loading factor down to our 
surface density of E = 10 M© pc“^ yields very similar re¬ 
sults. 


6 HIGHER AND LOWER SN RATES 

In all runs with lower SN rate (three times lower than the KS 
SN rate) the gas collapses into a very thin sheet. The energy 
input from the SNe is not able to support the disc against the 
gravitational attraction. Molecular gas forms very efficiently. 
As in the KS simulations, the peak and mixed SN models do 
not drive any outflow. Randomly placed SNe (SlO-lowSN- 


rand) are able to temporarily push gas out to 2 : ~ 700 pc, 
but at very low speeds. Only in the last few Myr of the 
simulation time a chimney establishes and a weak outflow is 
launched. 

In the case of high SN rates (three times as high as 
the KS SN rate) the evolution of the disc depends even 
more strongly on the driving mechanism. For random (SIO- 
highSN-rand) and mixed (SlO-highSN-mix) SN positions the 
kinetic pressure is high enough to completely evacuate the 
midplane region after ~ 20 — 30 Myr. As a factor of a few 
higher SN rate is not implausible for real systems, the tem¬ 
poral effects of starburst events and the resulting cavities in 
the disc midplane might not be unphysical, i.e. the existence 
of a temporally evacuated central region of the disc does not 
mark the model as unphysical. However, at this point the 
constant SN rate with positions in the evacuated region is 
not justihed any more. In simulation SlO-highSN-mix the 
gas remains in the region around zbO.bkpc with a balanced 
force from SN driving and gravity. Simulation SlO-highSN- 
rand does not show a turnover of the gas within the simu¬ 
lated time. In the case of SNe in density peaks (SlO-highSN- 
peak) the energy input is strong enough to efficiently mix 
the gas on small scales and thus prevents the formation of 
large density contrasts, so the disc has very smooth density 
distribution extending from —0.5 — 0.5 kpc. Similar to the 
run with the KS SN rate no outflow is launched. 

In Fig.[^we compare the velocity dispersions for peak, 
mixed, and random SN driving for low, KS, and high SN 
rates. The upper panel depicts crH 2 , the middle and lower 
panel present the observational estimates for Hi and Hq,. 
For crH 2 we notice an early trend with SN rate which be¬ 
comes weaker over time and is lost after ~ 30 Myr when the 
hrst clouds and voids in the disc have formed (see Fig. |^. 
The values span a range from 2 — 6kms“^. The evolution 
of am reveals an overall dependence on the SN rate but to¬ 
wards the end an even stronger dependence on the driving 
mechanism. SNe in density peaks lead to the lowest values 
followed by models with mixed and random SN position¬ 
ing. In (jHa the temporal scatter is very large but overall 
the velocity dispersions are determined to hrst order by the 
driving mechanism with the same ordering as for am- 


7 POTENTIAL CAVEATS 


7.1 Dynamical equilibrium 


Previous studies show that the formation of molecular gas 


in colliding flows of atomic gas takes about 10Myr (Clark 


et al. 2012). Two molecular clouds initially at rest with a 
mass oFTffi M© at a distance of 100 pc need about 40 Myr 
to merge due to their mutual gravitational attraction. This 
adds up to time-scales of around 50 Myr for the formation 
of a CMC. In a turbulent environment the formation and 
merging of clouds is likely to occur on shorter time-scales. 


which has also been suggested by observations (Tamburro 
et al.|20()8 Fukui et al.|2015 ) hnding time-scales as small as 


10" yr. 


Our simulated time-scale of 100 Myr is thus expected 
to cover the formation of molecular clouds, and indeed we 
observe this process in ah simulations (see Fig. where 
we show the total column density for simulation SIO-KS- 
rand over time). However, we see no dispersing or dissolving 
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Figure 19. Time evolution of the velocity dispersions for H 2 
(top) as well as the observational estimates for Hi (middle) and 
Hq; (top) for different SN rates and driving modes. A factor of 
three in the SN rate can double the velocity dispersion. Peak 
driving results in lower cr compared to mixed and random driving. 
At a later stage of the simulations the impact of the driving mode 
can be larger than a three times higher or lower SN rate. 


clouds. We conclude that SNe alone cannot destroy massive 
molecular clouds, even if they are all placed in density peaks, 
i.e. in the centre of molecular clouds. A dynamical equilib¬ 
rium in a sense of a constant number of clouds or a constant 
average cloud mass is not established in our simulation. 

This leads to the conclusion that either the assumption 
of a constant star formation and SN rate is not appropriate 
on the scales under consideration or that we neglect impor¬ 
tant physical processes in our model such as feedback. 


7.2 Assumption of a Constant SFR 

On large scales, i.e. when considering a signihcant fraction 
of the Galactic disc with a statistically relevant number of 
star-forming regions, the assumption of an averaged con¬ 


stant SN rate that scales directly with the surface den¬ 
sity is observationally well justihed. This is known as the 
Kennicutt-Schmidt (KS) relation ( Schmidt|1959 Kennicutt 


1998 


Bigiel et al. 2008). Locally, the KS relation breaks 


down, which is expected because individual star-forming 
complexes together with their temporal evolution need to 
be considered (see also Kruijssen Longmore|20T4 ). Time- 
dependent SN rates that are linked to the formation of 
molecular gas are likely to prevent the continuous increase 
of the average cloud mass by signihcantly increasing the lo¬ 
cal (spatial and temporal) SN feedback compared to the 
constant SN rate based on the KS relation. In simulations 
with very simplihed physics this problem might not occur at 
ah because of the missing stable dense gas structures that 
serve as star forming regions. Indeed, it has been shown that 
a simplihed thermodynamic model and no self-gravity does 
statistically result in a turbulent ISM with reasonable frac¬ 


tions of hot and cold gas (see, e.g. de Avillez & Breitschw- 


erdt 2004). However, no long-lived stable structures form 


in those simplihed setups. The small-scale transient features 
can still be treated and efficiently stirred with a globally ori¬ 
ented star formation and SN prescription to give reasonable 
values in terms of dynamics (e.g. velocity dispersion) and 
thermodynamics (e.g. volume-hlling fractions). 


7.3 Stellar feedback 


It is likely that additional stellar feedback processes pro¬ 
vide additional support against gravitational collapse, e.g. 
protestehar outflows, radiation feedback, and stellar winds. 
Early on during the formation process of massive stars pro- 
tostellar outflows can signihcantly stir the gas in molecular 
clouds ( Seifried et al.||20I2| Federrath et al.||20I4 Offner & 


Arce||2014 Peters et al.||20I4 ) providing mechanical energy 

against further collapse of molecular clouds. Radiation feed¬ 
back has long been suggested as an important heating source 
that can suppress gravitational collapse and expel gas from 
clouds ( Stromgr^|I939 Kahn||I954 Oort k, Spitzer||I955 ). 
Numerical studies support this paradigm and show efficient 


heating of the surroundings by OB stars (Dale et al. 2005 
Peters et al.|20I0l|Arthur et al.|20II||Walcli et al.|20I2 [Dale 

et al.|20I4[|Geen et al.|20I5| . Recent work by |Boneberg et al. 
(2015) shows that radiation feedback can even drive turbu¬ 


lence in molecular clouds. Although being generally less effi¬ 
cient than UV photoheating, direct radiation pressure is ex¬ 


pected to drive dynamical feedback (Krumholz & Matzner 
2009| [Murray, Quataert k, Thompson||20I0| [Krumholz k, 


Thompson 


2012 


[Sales et al.||2014 ). Stellar winds might 


equally well counteract gravitational attraction due to the 
formation of hot low-density bubbles, found by theoretical 


estimates ( 

Avedisova|I972 Gastor, McGray & Weaver|I975 

Weaver et al.||I977 

) as well as recent numerical simulations 

(Rogers & Pittard 

2013 Dale et al.||2014|. 


7.4 Global dynamics 

Our simulations also imply certain assumptions of the global 
dynamical properties of the galaxy under consideration. 
Most notably, we do not include the shear motions that are 
characteristics of the gas flow in typical disc galaxies. In or¬ 
der to determine the potential impact of galactic rotation 
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we evaluate the Rossby number, 


cr = 




(16) 


which is defined as the ratio of inertial to Coriolis force. 
Here, Ugas is the velocity of the gas, Q = Urot/(27rR) is the 
angular frequency with which the stratified box would rotate 
around the centre of the galaxy, Vrot the tangential velocity, 
R the distance from the galactic centre, and L is the size 
of the box in the rotation plane (dimension in x and y in 
our case). Thus small numbers indicate that the system is 
affected by the Coriolis force. For our box with L = 0.5 kpc 
and typical numb ers of Ugas ~ 10kms“^, Urot = 200kms“^ 
( Reid et al.|l999 ), R = 8kpc {Q ^ Sx 10“^^s“^) the Rossby 
number is cr ~ 5. The gas in our box reaches velocities a 
factor of a few higher and most of the structures we are con¬ 
centrating on are smaller than the size of the total box. Both 
corrections increase the Rossby number, so our simulations 
are not expected to be influenced by rotation. However, as 
shear might still be important for some parts of the dy¬ 
namics in the ISM, we are going to implement shearing box 
conditions for future projects. 


8 SUMMARY AND CONCLUSIONS 

With the SILCC simulations we aim at better understand¬ 
ing the dynamical evolution of the ISM and the life cycle of 
molecular clouds in star forming disc galaxies like the Milky 
Way. We simulate the evolution of a part of a galactic disc 
(a stratified box of size 0.5 x 0.5 x ±5 kpc with a gas sur¬ 
face density of E = 10 M© pc“^) on a ~ 100 Myr time-scale 
including an external stellar disc potential and self-gravity. 
Gas is heated by the background interstellar radiation and 
can cool based on a chemical network incorporating H^, H, 
H 2 , CO and C^, which also allows us to follow the formation 
of molecular gas. We include shielding of the gas and stellar 
feedback in the form of SNe. The SN rate is based on the 
Kennicutt-Schmidt (KS) relation assuming a universal ini¬ 
tial mass function. We cover the scatter in the KS relation by 
varying the SN rate by a factor 3 and 1/3. The positioning 
of the SNe is varied between random driving (random loca¬ 
tions), peak driving, where the SNe are placed in the density 
peaks, mixed driving (with equal random and peak contri¬ 
butions), and clustered driving, where a fraction of the SNe 
are spatially clustered following the available observational 
constraints. We also test the effect of magnetic fields with 
an initial field strength of H = 3 yG oriented in the plane of 
the disc. Some simulations form large molecular cloud com¬ 
plexes that become Jeans unstable and undergo collapse to 
form stars. At this time the calculation is stopped even if 
t < 100 Myr. 

Our main results can be summarized as follows: 

• The positioning of SN explosions plays a major role in 
determining the dynamics in the ISM. SNe placed at ran¬ 
dom positions lead to a higher gas velocity dispersion than 
SNe exploding in density peaks, where much of the input 
energy can be dissipated efficiently. The effect of SN posi¬ 
tioning might easily dominate over variations in the SN rate. 
Comparing velocity dispersions based on chemical composi¬ 
tion (H 2 ), observational estimates (Hi, Hq,) and tempera¬ 
ture regimes (8000 K; 3 x 10^ K, 300 K; 8000 K, < 300 K) 


shows good agreement with observations and previous ISM 
studies (ctri ~ 0 . 8 ( 7300-8000 k ~ 10 — 20kms cTRa ~ 
O-ScTgooo-sxios K ~ 20 — 30kms“^). We find almost con¬ 
stant ratios of the following velocity dispersions over time: 
^'Ki/cTRck 0.6, (7Ri/crR2 3 — 4, CTRa/cTgooO-SxlOS K ~ 
<7ri/(7300-8000K ~ 0.8, and (7R2/(7<300K ~ 0.4. 

• Randomly placed SNe (both individual SNe and ran¬ 
domly placed clusters of SNe) drive gaseous outflows after 
t ~ 10 — 30 Myr with mass loading factors of up to 10. In 
contrast, SNe placed in density peaks do not drive any no¬ 
ticeable outflow. Outflows are launched when the ISM is 
structured into filaments, clouds and voids and the SN rem¬ 
nants are able to expand into low-density gas. This allows 
for the formation of low-density chimneys and outflow chan¬ 
nels (~ 200pc wide), through which the gas is ejected with 
velocities of up to a few 100kms“^. The low-density, high- 
velocity gas {v > lOOkrns”^) is mostly ionized and drags 
denser atomic hydrogen with it. The bulk of the outflow¬ 
ing mass is dense (p ~ 10“^^ — 10“^^gcm“^) and slow 
(u ~ 20 — 40kms“^) but there is a high-velocity tail of 
up to u ~ 500kms“^ with p ~ 10“^® — 10“^^gcm“^, which 
significantly contributes to the total outflowing mass. 

• The outflows are predominantly composed of atomic 
hydrogen. If we model only clustered SNe the fraction of 
ionized hydrogen increases and the outflowing gas is com¬ 
posed of H and H^ in roughly equal proportions. Clus¬ 
tered SNe at random positions start driving outflows ear¬ 
lier than individual SNe at random positions. However, the 
overall outflow rates are very similar for both driving mecha¬ 
nisms. The primary driver for outflows is the vertical kinetic 
pressure gradient. In all simulations that generate outflows 
the kinetic pressure gradient is significantly larger than the 
thermal counterparts. Due to the finite simulation time of 
~ 100 Myr we are limited to studying the launching pro¬ 
cesses of outflows. Investigations of galactic fountains require 
longer time-scales and are beyond the scope of the paper. 

• The ISM evolves without reaching any kind of dynam¬ 
ical equilibrium over the entire simulation time of 100 Myr, 
i.e. dense filaments and small clumps merge to form GMCs. 
The clouds typically increase in size and mass. Due to their 
small volume filling fractions the massive dense clumps of 
molecular gas cannot be disrupted by randomly placed SNe 
and the peak SNe dissipate too much energy to efficiently 
disrupt the clouds. Depending on the positions of the SNe 
and the SN rate the formation of giant molecular clouds 
is abetted or retarded but a final coalescence of almost all 
the cold gas in one or a few giant clouds is inevitable. This 
is a potential caveat of our simulation setup as we do not 
attempt to include star formation self-consistently. We do 
follow the formation of molecular clouds but are unable to 
simulate their destruction. 

• Self-gravity is of major importance for developing a re¬ 
alistic structure of the disc. In the absence of self-gravity 
the amount of molecular gas becomes unrealistically low 
(/h 2 ^ 5% compared to /r 2 > 40% including self-gravity) 
and the velocity dispersion in the molecular gas drops to 
crR 2 ~ 2kms“^ ((7 r 2 ~ 5kms“^ including self-gravity). The 
densities and pressures in the midplane are at least an or¬ 
der of magnitude lower compared to simulations including 
self-gravity and compared to observed estimates. 

• Magnetic fields do not show a significant dynamical im¬ 
pact: the velocity dispersions, the pressures in the midplane 
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as well as the outflow properties are very similar to the cor¬ 
responding run without magnetic fields. At the end we find 
fewer but more massive clouds in the magnetic run, but this 
is likely to be influenced by the initial conditions for the 
magnetic field. At the beginning of the simulation the mag¬ 
netic tension efficiently reduces mixing and the magnetic 
pressure delays gravitational attraction which manifests in 
a delay of the formation of dense structures. Once molec¬ 
ular clouds have formed the simulations with and without 
magnetic fields behave very similarly. 

We find the best match to the observational data with 
randomly placed individual or clustered supernovae. Our 
calculations demonstrate that a magnetised, self-gravitating 
disc-like sheet of gas driven at a rate consistent with the 
global KS relation evolves on relatively short time-scales of 
~ 50Myr into a ’realistic’ multi-phase ISM. Its morphologi¬ 
cal and kinematic structure as well as its phase fractions and 
outflow properties agree well with the observational data. 
Future investigations focusing on a more self-consistent star 
formation treatment and longer time-scales will be required 
to establish a potential ’self-regulated’ or ’equilibrium’ state 
and estimate the global impact on galaxy evolution. 
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